/*
 * Calculates calibration, quadratic fit.
 *
 * ge_calib_e_quad spams calib_single with guesses for where a specific peak
 * is. calib_single guesses, looks for peaks and does LSQF.
 */

#include <iomanip>

std::ofstream file;

void ge_calib_e_quad()
{
  file.open("out.dat");
  calib_single(EGamCl0Cr0, 11422);
  calib_single(EGamCl0Cr0Se0, 11423);
  calib_single(EGamCl0Cr0Se1, 11441);
  calib_single(EGamCl0Cr0Se2, 11430);
  calib_single(EGamCl0Cr0Se3, 11517);
  calib_single(EGamCl0Cr0Se4, 11426);
  calib_single(EGamCl0Cr0Se5, 11412);
  calib_single(EGamCl0Cr1, 11430);
  calib_single(EGamCl0Cr1Se0, 11440);
  calib_single(EGamCl0Cr1Se1, 11427);
  calib_single(EGamCl0Cr1Se2, 11433);
  calib_single(EGamCl0Cr1Se3, 11436);
  calib_single(EGamCl0Cr1Se4, 11430);
  calib_single(EGamCl0Cr1Se5, 11420);
  calib_single(EGamCl0Cr2, 11420);
  calib_single(EGamCl0Cr2Se0, 11426);
  calib_single(EGamCl0Cr2Se1, 11424);
  calib_single(EGamCl0Cr2Se2, 11428);
  calib_single(EGamCl0Cr2Se3, 11434);
  calib_single(EGamCl0Cr2Se4, 11427);
  calib_single(EGamCl0Cr2Se5, 11426);
  calib_single(EGamCl1Cr0, 11419);
  calib_single(EGamCl1Cr0Se0, 11418);
  calib_single(EGamCl1Cr0Se1, 11419);
  calib_single(EGamCl1Cr0Se2, 11410);
  calib_single(EGamCl1Cr0Se3, 11439);
  calib_single(EGamCl1Cr0Se4, 11423);
  calib_single(EGamCl1Cr0Se5, 11426);
  calib_single(EGamCl1Cr1, 11418);
  calib_single(EGamCl1Cr1Se0, 11422);
  calib_single(EGamCl1Cr1Se1, 11418);
  calib_single(EGamCl1Cr1Se2, 11418);
  calib_single(EGamCl1Cr1Se3, 11419);
  calib_single(EGamCl1Cr1Se4, 11416);
  calib_single(EGamCl1Cr1Se5, 11428);
  calib_single(EGamCl1Cr2, 11422);
  calib_single(EGamCl1Cr2Se0, 11426);
  calib_single(EGamCl1Cr2Se1, 11426);
  calib_single(EGamCl1Cr2Se2, 11421);
  calib_single(EGamCl1Cr2Se3, 11438);
  calib_single(EGamCl1Cr2Se4, 11406);
  calib_single(EGamCl1Cr2Se5, 11429);
  calib_single(EGamCl2Cr0, 11426);
  calib_single(EGamCl2Cr0Se0, 11451);
  calib_single(EGamCl2Cr0Se1, 11432);
  calib_single(EGamCl2Cr0Se2, 11437);
  calib_single(EGamCl2Cr0Se3, 11425);
  calib_single(EGamCl2Cr0Se4, 11432);
  calib_single(EGamCl2Cr0Se5, 11425);
  calib_single(EGamCl2Cr1, 11421);
  calib_single(EGamCl2Cr1Se0, 11417);
  calib_single(EGamCl2Cr1Se1, 11424);
  calib_single(EGamCl2Cr1Se2, 11417);
  calib_single(EGamCl2Cr1Se3, 11402);
  calib_single(EGamCl2Cr1Se4, 11411);
  calib_single(EGamCl2Cr1Se5, 11415);
  calib_single(EGamCl2Cr2, 11420);
  calib_single(EGamCl2Cr2Se0, 11413);
  calib_single(EGamCl2Cr2Se1, 11419);
  calib_single(EGamCl2Cr2Se2, 11487);
  calib_single(EGamCl2Cr2Se3, 11417);
  calib_single(EGamCl2Cr2Se4, 11416);
  calib_single(EGamCl2Cr2Se5, 11422);
  calib_single(EGamCl3Cr0, 11427);
  calib_single(EGamCl3Cr0Se0, 11440);
  calib_single(EGamCl3Cr0Se1, 11414);
  calib_single(EGamCl3Cr0Se2, 11420);
  calib_single(EGamCl3Cr0Se3, 11425);
  calib_single(EGamCl3Cr0Se4, 0);
  calib_single(EGamCl3Cr0Se5, 11424);
  calib_single(EGamCl3Cr1, 11422);
  calib_single(EGamCl3Cr1Se0, 11431);
  calib_single(EGamCl3Cr1Se1, 11427);
  calib_single(EGamCl3Cr1Se2, 11410);
  calib_single(EGamCl3Cr1Se3, 11417);
  calib_single(EGamCl3Cr1Se4, 0);
  calib_single(EGamCl3Cr1Se5, 11424);
  calib_single(EGamCl3Cr2, 11424);
  calib_single(EGamCl3Cr2Se0, 0);
  calib_single(EGamCl3Cr2Se1, 11420);
  calib_single(EGamCl3Cr2Se2, 11428);
  calib_single(EGamCl3Cr2Se3, 11416);
  calib_single(EGamCl3Cr2Se4, 0);
  calib_single(EGamCl3Cr2Se5, 11443);
  calib_single(EGamCl4Cr0, 11426);
  calib_single(EGamCl4Cr0Se0, 11432);
  calib_single(EGamCl4Cr0Se1, 11420);
  calib_single(EGamCl4Cr0Se2, 11430);
  calib_single(EGamCl4Cr0Se3, 11428);
  calib_single(EGamCl4Cr0Se4, 11429);
  calib_single(EGamCl4Cr0Se5, 11454);
  calib_single(EGamCl4Cr1, 11419);
  calib_single(EGamCl4Cr1Se0, 11417);
  calib_single(EGamCl4Cr1Se1, 11423);
  calib_single(EGamCl4Cr1Se2, 11426);
  calib_single(EGamCl4Cr1Se3, 11419);
  calib_single(EGamCl4Cr1Se4, 11415);
  calib_single(EGamCl4Cr1Se5, 11427);
  calib_single(EGamCl4Cr2, 11439);
  calib_single(EGamCl4Cr2Se0, 11444);
  calib_single(EGamCl4Cr2Se1, 11439);
  calib_single(EGamCl4Cr2Se2, 11449);
  calib_single(EGamCl4Cr2Se3, 11440);
  calib_single(EGamCl4Cr2Se4, 11440);
  calib_single(EGamCl4Cr2Se5, 11462);
  calib_single(EGamCl5Cr0, 11420);
  calib_single(EGamCl5Cr0Se0, 11422);
  calib_single(EGamCl5Cr0Se1, 0);
  calib_single(EGamCl5Cr0Se2, 11421);
  calib_single(EGamCl5Cr0Se3, 11420);
  calib_single(EGamCl5Cr0Se4, 11427);
  calib_single(EGamCl5Cr0Se5, 0);
  calib_single(EGamCl5Cr1, 11424);
  calib_single(EGamCl5Cr1Se0, 11438);
  calib_single(EGamCl5Cr1Se1, 11423);
  calib_single(EGamCl5Cr1Se2, 11433);
  calib_single(EGamCl5Cr1Se3, 11435);
  calib_single(EGamCl5Cr1Se4, 0);
  calib_single(EGamCl5Cr1Se5, 11414);
  calib_single(EGamCl5Cr2, 11432);
  calib_single(EGamCl5Cr2Se0, 11429);
  calib_single(EGamCl5Cr2Se1, 11435);
  calib_single(EGamCl5Cr2Se2, 11439);
  calib_single(EGamCl5Cr2Se3, 0);
  calib_single(EGamCl5Cr2Se4, 11444);
  calib_single(EGamCl5Cr2Se5, 11452);
  calib_single(EGamCl6Cr0, 11409);
  calib_single(EGamCl6Cr0Se0, 11417);
  calib_single(EGamCl6Cr0Se1, 11416);
  calib_single(EGamCl6Cr0Se2, 11415);
  calib_single(EGamCl6Cr0Se3, 11411);
  calib_single(EGamCl6Cr0Se4, 11409);
  calib_single(EGamCl6Cr0Se5, 11415);
  calib_single(EGamCl6Cr1, 11418);
  calib_single(EGamCl6Cr1Se0, 11416);
  calib_single(EGamCl6Cr1Se1, 11413);
  calib_single(EGamCl6Cr1Se2, 11415);
  calib_single(EGamCl6Cr1Se3, 11416);
  calib_single(EGamCl6Cr1Se4, 11409);
  calib_single(EGamCl6Cr1Se5, 11411);
  calib_single(EGamCl6Cr2, 11412);
  calib_single(EGamCl6Cr2Se0, 0);
  calib_single(EGamCl6Cr2Se1, 11395);
  calib_single(EGamCl6Cr2Se2, 11413);
  calib_single(EGamCl6Cr2Se3, 11418);
  calib_single(EGamCl6Cr2Se4, 11399);
  calib_single(EGamCl6Cr2Se5, 11413);
  calib_single(EGamCl7Cr0, 11414);
  calib_single(EGamCl7Cr0Se0, 11433);
  calib_single(EGamCl7Cr0Se1, 11416);
  calib_single(EGamCl7Cr0Se2, 11428);
  calib_single(EGamCl7Cr0Se3, 11422);
  calib_single(EGamCl7Cr0Se4, 11421);
  calib_single(EGamCl7Cr0Se5, 11428);
  calib_single(EGamCl7Cr1, 11423);
  calib_single(EGamCl7Cr1Se0, 11436);
  calib_single(EGamCl7Cr1Se1, 11426);
  calib_single(EGamCl7Cr1Se2, 11423);
  calib_single(EGamCl7Cr1Se3, 11424);
  calib_single(EGamCl7Cr1Se4, 11423);
  calib_single(EGamCl7Cr1Se5, 11434);
  calib_single(EGamCl7Cr2, 11428);
  calib_single(EGamCl7Cr2Se0, 11134);
  calib_single(EGamCl7Cr2Se1, 11423);
  calib_single(EGamCl7Cr2Se2, 11396);
  calib_single(EGamCl7Cr2Se3, 11429);
  calib_single(EGamCl7Cr2Se4, 11418);
  calib_single(EGamCl7Cr2Se5, 11431);
  file.close();
}

static double subdet(double m[3][3], int row, int col)
{
/*  int row1, col1, row2, col2;
  switch (row) {
    case 0: row1 = 1; row2 = 2; break;
    case 1: row1 = 0; row2 = 2; break;
    case 2: row1 = 0; row2 = 1; break;
  }
  switch (col) {
    case 0: col1 = 1; col2 = 2; break;
    case 1: col1 = 0; col2 = 2; break;
    case 2: col1 = 0; col2 = 1; break;
  }*/
  int row1 = (2 - row) / 2;
  int row2 = (5 - row) / 2;
  int col1 = (2 - col) / 2;
  int col2 = (5 - col) / 2;
  return m[row1][col1] * m[row2][col2] - m[row1][col2] * m[row2][col1];
}

static void calib_single(TH1F const *hist, int guess_21Ne)
{
  static int ge_id = 0;

  int const kPeakNum = 8;
  double const kPeakE[kPeakNum] = {122, 245, 344, 779, 965, 1087, 1113, 1409};
  int const kPeakWidth = 40;

  double guess_k = guess_21Ne / kPeakE[7];
  double peak_e[kPeakNum];

  // Guess peaks.
  for (int peak_idx = 0; peak_idx < kPeakNum; ++peak_idx) {
    int const kGuess = kPeakE[peak_idx] * guess_k;
    int peak_max = 0;
    int peak_max_i = 0;
    for (int i = kGuess - kPeakWidth/2; i < kGuess + kPeakWidth/2; ++i)
      if (hist->GetBinContent(i) > peak_max) {
        peak_max = hist->GetBinContent(i);
        peak_max_i = i;
      }
    double sum_fx = 0.;
    double sum_f = 0.;
    for (int i = peak_max_i - kPeakWidth/2; i < peak_max_i + kPeakWidth/2; ++i)
    {
      sum_fx += i * hist->GetBinContent(i);
      sum_f += hist->GetBinContent(i);
    }
    peak_e[peak_idx] = sum_fx / sum_f;
  }

  double X[kPeakNum][3];
  for (int row = 0; kPeakNum > row; ++row) {
    X[row][0] = peak_e[row] * peak_e[row];
    X[row][1] = peak_e[row];
    X[row][2] = 1.0;
  }
  double X_T_X[3][3];
  for (int row = 0; 3 > row; ++row) {
    for (int col = 0; 3 > col; ++col) {
      double sum = 0.0;
      for (int i = 0; kPeakNum > i; ++i) {
        sum += X[i][row] * X[i][col];
      }
      X_T_X[row][col] = sum;
    }
  }
  double idet = 1.0 / (
      + X_T_X[0][0] * X_T_X[1][1] * X_T_X[2][2]
      + X_T_X[0][1] * X_T_X[1][2] * X_T_X[2][0]
      + X_T_X[0][2] * X_T_X[1][0] * X_T_X[2][1]
      - X_T_X[0][2] * X_T_X[1][1] * X_T_X[2][0]
      - X_T_X[0][1] * X_T_X[1][0] * X_T_X[2][2]
      - X_T_X[0][0] * X_T_X[1][2] * X_T_X[2][1]
      );
  double X_T_X_1[3][3];
  for (int row = 0; 3 > row; ++row) {
    for (int col = 0; 3 > col; ++col) {
      X_T_X_1[row][col] = idet * subdet(X_T_X, row, col) *
          (0 == (1 & (col + row * 3)) ? 1 : -1);
    }
  }

  double X_T_y[3];
  for (int row = 0; 3 > row; ++row) {
    double sum = 0.0;
    for (int i = 0; kPeakNum > i; ++i) {
      sum += X[i][row] * kPeakE[i];
    }
    X_T_y[row] = sum;
  }
  double B[3];
  for (int row = 0; 3 > row; ++row) {
    double sum = 0.0;
    for (int i = 0; 3 > i; ++i) {
      sum += X_T_X_1[row][i] * X_T_y[i];
    }
    B[row] = sum;
  }

  int clu_id = ge_id / (7 * 3);
  int cry_id = (ge_id / 7) % 3;
  int seg_id = ge_id % 7;
  int dgf_id = 6 * clu_id + 2 * cry_id + (3 > seg_id ? 0 : 1);
  int sign_id = 3 > seg_id ? seg_id : seg_id - 3;
  std::ostringstream oss;
  oss << "dgf_" << dgf_id << '_' << sign_id << '.';
  file << oss.str() << "Gain: " << B[1] << '\n';
  file << oss.str() << "Offset: " << B[2] << '\n';

  ++ge_id;
}
